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1  Introduction 

Under  this  grant  we  have  attemted  to  simulate  supersonic  reactive  flows  with  high  order 
accuracy  methods.  The  scientific  reason  for  using  high  order  methods  is  that  simulations 
of  supersonic  reactive  flows  require  long  time  integrations  and  the  resolution  of  fine  scales 
of  the  flow.  It  is  well  known  that  high  order  accuracy  methods  are  mandatory  in  satisfying 
these  requirements. 

Our  main  codes  are  spectral,  i.e.  based  on  expansions  in  global  polynomials.  The  two 

issues  that,  we  FAve_  adm:oHsed.lmAhA 

flexibility.  To  adapt  spectral  methods  to  discontinuities  we  used  low  pass  filters  to  stabilize 
the  scheme  and  a  postprocessing  of  the  solution  to  recover  the  accuracy.  We  have  resolved 
the  Gibbs  phenomenon  and  showed  that  a  proper  postprocessing  can  recover  spectral 
accuracy  in  smooth  regions  of  the  flow.  The  theory  is  not  a  constructive  one  in  the  sense 
that  it  does  not  discuss  optimal  methods  for  postprocessing. 

A  great  deal  of  work  has  been  done  in  the  last  decade  in  adaptive  spectral  methods  to 
complex  domains.  The  idea  here  is  to  use  multidomain  techniques  where  spectral  methods 
are  used  locally  at  every  subdomain.  Good  sets  of  collocation  points  has  been  found  for 
triangles  and  tetrahydra,  that  enable  the  construction  of  interpolation  polynomials  in  those 
subdomains.  The  main  difficulty  is  the  imposition  of  interface  boundary  conditions.  We 
have  started  to  study  an  alternative  of  the  DG  method.  We  look  for  a  discontinuous 
Collocation  method,  which  4s  more  natural  in  the  framework  of  spectral  methods. 

There  are  uncertainties  associate  with  any  aspect  of  computations  of  reactive  flows. 
We  started  to  study  the  effects  of  uncertainties  in  initial  conditions  on  the  steady  state 


supersonic  flows  in  the  double-throated  nozzle  using  the  polynomial  chaos  expansions.  The 
results  obtained  are  summarized  in  Section  4. 

In  the  applications  part  we  considered  two  major  applications:  The  first  is  simulations 
of  Shock-Induced  Turbulence  Mixing  and  the  second  is  the  computational  study  of  flame- 
holders  in  SCRAMJETs,  mainly  recessed  cavity  flameholders.  For  the  first  subject  the 
principal  research  objective  has  been  to  develop  new  state-of-the-art  high-order  accurate 
numerical  methods  for  the  multi-dimensional  numerical  simulations  of  the  fully-nonlinear 
evolution  of  hydrodynamic  instabilities  and  late-time  turbulent  mixing  generated  by  single- 
and  multi-mode  Richtmyer-Meshkov,  Rayleigh- Taylor,  and  Kelvin-Helmholtz  instabilities. 
This  will  be  discussed  in  Section  6. 

In  the  second  application  we  considered  supersonic  combustion  problems  in  recessed 
cavities  in  order  to  establish  the  efficacy  of  recessed  cavity  flame-holders.  Recessed  cavities 
provide  a  high  temperature,  low  speed  recirculating  region  that  can  support  the  produc¬ 
tion  of  radicals  created  during  chemical  reactions.  This  stable  and  efficient  flame-holding 
performance  by  the  cavity  is  achieved  by  generating  a  recirculation  region  inside  the  cav¬ 
ity  where  a  hot  pool  of  radicals  forms  resulting  in  reducing  the  induction  time  and  thus 
obtaining  the  auto-ignition.  Experiments  have  shown  that  such  efficiency  depends  on  the 
geometry  of  the  cavity  such  as  the  degree  of  the  slantness  of  the  aft  wall  and  the  length  to 
depth  ratio  of  cavity.  This  will  be  discussed  in  Section  7. 


2  Postprocessing  Techniques 


Sufficient  conditions  fox  the  removal  of  the  Gibus'  phenomenon  were  gi  veil  hi  [38] .  Consider 
a  function  f(x)  €  L2[— 1, 1]  and  assume  that  there  is  a  subinterval  [a,b]  C  [—1, 1]  in  which 
f(x )  is  analytic.  (For  convenience  we  define  the  local  variable,  £  =  —  1  +  such  that 

if  a  <  x  <  b  then  —  1  <  £  <  1.)  Let  the  family  {^(a;)},  be  orthonormal  under  a  scalar 
product  ( • ,  • ),  and  denote  the  finite  expansion  of  f(x)  in  this  basis  by  /w(:r), 


N 


In(x)  = 

k—0 


Let  the  family  {$£(£)}  be  Gibbs  complementary  to  the  family  {4t(x)}  (see  [38]  for  its 
exact  definition),  then  the  postprocessed  reconstruction  given  by 

g«(x)  =  Z  >a  #?(«*)) 

1=0 

converges  exponentially  to  f(x),  i.e. 


max  \f{x)  -  gN(x)\  < 

a<x<b 


„-qN 


q  >  0. 


9 


In  a  series  of  papers  it  has  been  shown  that  the  Gegenbauer  polynomials 

*1(0  =  4rC(A(f) 

which  are  orthonormal  under  the  inner  product  <  •,  •  >A  defined  by 

<f,9>  a  =  J_jl-a?)x->f{£)g(t)d£ 

are  Gibbs  complementary  to  all  commonly  used  spectral  approximations. 

The  Gegenbauer  method  is  not  robust,  it  is  sensitive  to  roundoff  errors  and  to  the  choice 
of  the  parameters  A  and  m.  A  different  implementation  of  the  Gegenbauer  postprocessing 
method  has  been  suggested  recently  by  Jung  and  Shizgal  [53].  To  explain  the  differences 
between  the  direct  Gegenbauer  method  and  the  inverse  Gegenbauer  method  suggested  in 
[53],  consider  the  case  of  the  Fourier  expansion  of  a  nonperiodic  problem.  The  Fourier 
approximation  /^(x)  of  f(x) 

Mx)  =  £  Aeite, 

k=—N 

where  fk  =  (/(x),elfc7ri),  and  we  construct 

m 

fS(x)  =  £j,cA(*), 

1=0 

where  gi  =<  //v,  Cx(x)  >>.  In  the  Inverse  method  we  use  t.he  relationship 


h  =  (Jtf 


and  invert  to  find  <&.  Thus  if  we  define  the  matrix 

Wkl  =  {C?(x),eik™)F  =  p  i  1  -  x2)x~1W^{x)eiknxdx, 
and  fk  =  (/,  elk7rx),  then 

m 

£  =  ft- 

1=0 

The  method  seems  to  be  less  sensitive  to  roundoff  errors  or  to  the  choice  of  parameters. 
In  particular  if  the  original  function  is  a  polynomial,  the  inverse  method  is  exact. 

Another  approach  for  the  post  processing  issue  involves  rational  functions.  A  rapidly 
converging  approximation  to  finite  Fourier  series  has  been  studied  by  defining  a  rational 
function  which  denominator  and  numerator  are  represented  in  finite  Fourier  sum.  Pade 
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rational  approximations  have  been  considered  as  one  of  the  popular  computational  methods 
of  representing  functions,  especially  rapidly  converging  functions  since  it  was  proposed  by 
H.  Pade  in  1892  [8,  49].  They  are  generally  more  efficient  than  polynomial  approximations 
with  a  reduced  number  of  operations  for  the  same  accuracy  [23,  25,  26]. 

The  methods  discussed  before  require  the  knowledge  of  the  position  of  discontinuity, 
however  with  no  knowledge  of  the  singularity,  Pade  reconstruction  recovers  back  a  non- 
oscillatory  solution  successfully  with  a  reduced  overshoot  at  the  singularity.  It  is  due  to 
the  fact  that  the  possible  existence  of  poles  of  some  order  for  the  denominator  of  Pade 
approximant  allows  to  give  a  better  approximation  to  the  functions  exhibiting  singular 
behaviors  such  as  large  gradient  and  discontinuity. 

Geer  and  his  coworkers  [26],  suggested  a  way  of  implementing  the  rational  trigonometric 
approximations  for  even  or  odd  27r-periodic  piecewise  smooth  functions  and  also  a  way  of 
applying  this  method  to  the  solution  of  the  initial  boundary  value  problem  heat  equation. 
In  their  work,  Fourier-Pade  approximants  are  defined  in  a  nonlinear  way  such  that  the 
relation  between  the  coefficients  of  the  rational  approximations  and  the  Fourier  coefficients 
involves  a  necessary  procedure  of  calculating  the  integration  of  rational  functions,  which 
makes  the  numerical  scheme  relatively  complicated.  In  [25],  Fourier  expansion  is  treated  as 
a  Laurent  expansion,  and  using  Fourier-Pade  rational  approach  the  spectral  convergence  is 
obtained  up  to  the  discontinuity  by  subtracting  off  the  jump  from  the  Fourier  data,  which 
requires  the  advance  knowledge  not  only  of  the  position  of  the  singularity,  but  also  the 
magnitude  of  the  jump. 

In  our  work,  we  have  designed  two  Fourier-Pade  methods  considering  the  general  case  of 
piecewise  analytic  functions,  with  no,  .adyauce jtmrwIedge^pLthe  singularity.:  Simple  ways  of. 
implementing  Fourier-Pade  Galerkin  and  Fourier-Pade  collocation  methods  are  developed 
and  applied  to  simulate  the  solutions  of  nonlinear  partial  differential  equations.  For  the 
hyperbolic  partial  differential  equations  such  as  Burgers’  equation,  an  initially  smooth 
function  can  evolve  into  shock  in  time  for  inviscid  case  and  large  gradient  for  viscous  case. 
Therefore  the  standard  spectral  simulations  exhibit  the  Gibbs  phenomenon  that  degrade 
the  accuracy  of  the  numerical  solutions  in  time.  From  accurate  Fourier  data  computed  by 
Fourier  method,  we  applied  the  Fourier-Pade  reconstruction  as  a  post- processing.  After  the 
post-processing,  the  computational  results  show  successful  reduction  of  the  non-physical 
oscillations  in  the  standard  spectral  solutions  of  the  one  dimensional  inviscid  Burgers’ 
equation  and  the  two  dimensional  inviscid  Boussinesq  convection  flows.  Especially,  the 
numerical  results  for  the  Boussinesq  convection  equation  are  demonstrated  in  figures  20 
and  23.  Further  study  has  to  be  established  to  find  the  optimal  relation  between  the  degrees 
of  the  polynomial  of  the  Pade  approximants  and  the  number  of  the  Fourier  coefficients. 
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3  Multi-Domain  Methods 


To  give  spectral  methods  flexibility  to  handle  complex  geometries  we  advocate  the  use  of 
multidomain  methods  with  penalty  type  interface  conditions,  i.e.  the  boundary  conditions 
are  imposed  only  in  a  weak  form. 

We  consider  two  interface  conditions,  i.e. 


1.  The  averaging  method,  in  which  the  interface  conditions  are  obtained  by  averaging 
the  state  vectors  of  the  two  adjacent  domains,  and 

2.  The  Penalty  method  in  conservative  form  in  which  the  interface  conditions  are  sat¬ 
isfied  only  in  a  weak  form,  leaving  the  approximations  not  necessarily  continuous  at 
the  interfaces. 


In  the  following  we  will  give  the  penalty  interface  conditions  for  the  Euler  and  Navier- 
Stokes  equations  and  also  show  that  the  averaging  method  is  a  subset  of  the  penalty 
method. 

Consider  the  inviscid  part  only,  in  the  ^-direction  in  the  interval  —  1  <  x  <  1,  i.e., 


dq  a F 

dt  dx 


=  0. 


(1) 


For  simplicity,  assume  that  we  have  two  domains  in  this  interval  with  the  interface  at 
x  =  0,  qIN{x ,  t )  denotes  the  numerical  solution  in  the  left  domain  x  <  0  and  q^ix,  t)  in  the 
right  domain  x  >  0.  Note  that  the  numerical  solution  is  composed  of  two  polynomials  of 

different  orders.  Tfeedbegeadre  •  fey  . •  -• 


dt  dx 


dt  dx 


M))  + 

riQAr(x)[/+(^(0,t))  -  /+(9m(M))]  + 
72Qw(ic)[/“(9w(0,<))  -  /~(?m(M))], 
B{qM(l,t))  + 

r3QM{x)[f+{qM{0,t ))  -  /+(^(0,f))]  + 

T4QM{x)[f~{q^{0,t))  -  f~{qlN{0,t))] 


(2) 


where  B  is  a  boundary  operator  at  the  end  points,  i.e.,  x  =  ±1  and  PN  and  Ijfj  are  the 
Legendre  interpolation  operators  for  the  left  and  right  domains  respectively.  .  The  positive 
and  negative  fluxes  f+  and  f~  are  defined  by 

f±  =  J  SA±S~1dq, 


F, 


(3) 


with 


(4) 


SAS~l. 


The  Jacobian  matrix  A  is  assumed  to  be  symmetric.  A+  and  A-  are  the  diagonal  matrices 
composed  of  positive  and  negative  eigenvalues  of  A  respectively.  Qn(x)  and  Qm(%)  are 
polynomials  of  orders  N  and  M  respectively  such  that  they  are  zero  at  all  the  collocation 
points  except  the  interface  points  x  =  0  (for  example  Qn(x)  =  — 0  <  x  <  1 
where  TN(x)  is  the  Chebyshev  polynomial  of  degree  N).  The  penalty  parameters  n,  72,73 
and  74  are  all  constants.  Since  we  are  interested  only  in  the  interface  conditions,  we 
ignore  the  boundary  operator  B  at  x  —  ±1.  Define  the  discrete  scalar  product  {p,g)N  = 
T,iLoPT(^i)<l(^i)u;i-  Ui  is  the  weight  in  the  Gauss-Lobatto-Legendre  quadrature  formula. 
With  the  discrete  product,  the  energy  Eft)  is  defined  by  Eft )  =  (qi/(x,t),qIN(x,t))ff  + 
(q^fx,  t),  q^fx,  t))M-  The  stability  conditions  of  penalty  parameters  are  given  by  the 
following  theorem: 


Theorem  1  The  energy  is  bounded  by  the  initial  energy  of  the  system  if  the  following 
conditions  are  satisfied  ; 

2 u)*nti  <  1,  2ujjnt2  >  1,  2 oj^t3  <  -1,  2 >  -1, 

UnT 1  -  w^r3  =  1,  a ;^r2  -  =  1.  (5) 


The  penalty  method  in  the  case  of  the  2-D  Euler  equation  is  given  by 


difi 

dt  dx  dy 


=  n ,3Q{x,y)[f+{qN)  -  f+(qM-)]  + 

72,4 Q(x,y)[f~(qN)  -  f~{qM-)]i 


where  qu-  is  the  state  vector  of  the  adjacent  domain  at  the  interface  of  degree  M,  ri)3(r2)4) 
denotes  rx (t2)  and  r3(r4)  respectively.  rx  and  r2  (r3  and  r4)  are  the  penalty  parameters 
for  the  right(left)  in  x-direction  and  top(bottom)  in  y-direction  respectively.  Q(x,  y)  is  a 
polynomial  which  vanishes  at  all  of  interior  points  of  the  domain  and  is  equal  to  1  at  the 
four  interfaces.  Note  that  the  boundary  operator  B  does  not  appear  in  the  scheme.  Let  A 
be  the  linearized  Jacobian  matrix  (around  a  state  vector  qf)  of  two  inviscid  fluxes 


(dF_  dG\ 
\dq ’  dq) 


(7) 


where  n  =  ( nx,ny )  is  the  unit  outward  normal  vector.  Since  the  matrix  A  is  symmetric, 
there  exists  S  such  that 


A  =  SAS~\ 


(8) 


where  A  is  a  diagonal  matrix  composed  of  eigenvalues  of  A.  Then  A  =  A+  +  A  and 
A±  =  SA±S~1.  A1*  is  defined  as  in  previous  section.  Splitting  A  yields 

/±  =  A±q0, .  (9) 

where  f±  is  obtained  from  the  linearized  state. 

Remark  1  Since  ft  =  (nx,  ny)  is  taken  to  be  outward  normal  vector,  the  stability  condition 
(6)  is  now  given  by 


2 caffTi  <  1,  2u!nT2  >  1,  2 w^r3  <  1,  2a$T4  >  1, 

+  wmt4  =  1,  oonT2  +  =  1.  (10) 


When  dealing  with  the  Navier  Stokes  equation,  we  keep  the  penalty  form  for  the  Euler 
fluxes  and  add  a  penalty  term  for  the  viscous  fluxes.  The  stability  of  this  procedure  stems 
from  the  fact  that  the  Jacobian  matrices  for  the  full  reactive  Navier-Stokes  equation  can 
be  symmetrized  by  the  same  similarity  transformation  (see  Appendix_B  in  [17]).  Thus  we 
get  the  system: 

dqM  dlfjF  dING 

dt  dx  dy 


I- ' '  :: . . (ii) 


dINFu  dING„ 

T\,sQ{x,y)[f+(qN )  -  f+{qM -)]  + 
T2,4Q{x,y)[f~{qN )  -  f~(QM -)]  + 
nfiQ(x,y)[Kv  qN-A„-  qM_]  + 


Here  /*  are  same  as  defined  in  the  previous  section  and  the  Jacobian  matrix  vector  A„  is 
given  by 


A„  = 


(dF„ 

\  dqx  Hx 


90  5 


(12) 


and 


Q  =(Q,Q),  dq  =  (qx,qy), 


(13) 


where  again  q_  and  3q_  denote  the  adjacent  domains  state  vectors  and  their  derivatives. 
Note  that  the  penalty  terms  A„  •  dq  does  not  appear  in  [9,  47,  48].  The  penalty  parameters 
t5) 7  and  76,8  are  defined  in  the  same  way  as  in  the  previous  section.  To  seek  stable  penalty 
parameters  we  split  the  inviscid  and  viscous  fluxes  and  keep  the  stability  conditions  of 
ti,2,3,4  for  the  inviscid  flux  as  in  Theorem  1.  The  stability  conditions  of  r5)7  and  r6)8  are 
given  in  the  following  Theorem  : 
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Theorem  2  The  penalty  method  for  the  Navier-Stokes  equations  (12)  is  stable  if  the 
penalty  parameters  tj  ,  j  =  1...4  are  as  in  Theorem  1  and  the  rest  satisfy: 

WNr6  <  0, 

—  OJmTs  =  0, 

1  +  unt5  -  umt7  =  0, 

—  +  — )  w2mt7  -  2 r7  +  4ujnTs  +  —  <0  .  (14) 

%  OJN  /  % 

in  addition  to  the  conditions  specified  for  t\,T2,  73  and  in  Theorem  1. 

Note  that  these  conditions  are  given  independently  of  the  local  flow  properties.  And 
moreover,  the  penalty  parameters  of  each  domain  are  constrained  by  its  adjacent  domain. 

Remark  2  For  n  to  be  outward  normal  vector  the  condition  (16)  is  now  given  by 

<  0,  WnTq  +  CJA/Tg  =0,  1  +  +  UJmT 7  =  0, 

- 1 - ^  u^r2  +  27>  +  H - <  0  (15) 

LOm  OJn  /  OJm 

with  the  conditions  (11) 

The  averaging  method  for  the  N-S  equations  can  be  presented  as 
dq  OF  dG  dFv  dGu 

-  -j-  -  — U  -  -  — j—  -  4- 

dt  dx  dy  dx  dy 

T2AQ{x,y)[f'~ {q)  -f'~(q-)]  + 
n,7 Q{x,  y)[Au  •  d?q  -  A„  •  <92q_]  + 
t6>8Q(x,  y)[Au  ■  d<\-  Au-  dq_], 

where  d2q  is  the  second  derivative  of  q  in  either  x  or  y  direction. 

Theorem  3  If  Ti  =  r3  =  7*2  =  r4  =  r5  =  r7  =  and  r6  =  — r8  =  —  then  the 

approximation  is  continuous  at  the  interface  and  the  scheme  (17)  is  stable. 

To  ensure  the  stability  of  the  scheme  at  some  particular  collocation  points  where  the 
solution  become  singular  and  unstable,  we  use  the  averaging  method  adaptively  at  selective 
grid  points.  In  particular,  we  switched  from  the  penalty  method  to  the  averaging  when  the 
following  criteria  was  satisfied: 


where  Cave  is  a  non-negative  constant.  Note  that  Cave  =  0  leads  to  the  averaging  method, 
whereas  a  large  Cave  results  in  the  penalty  method.  For  the  value  of  Cave  used  in  this 
paper,  we  found  out  that  there  were  very  few  points  in  which  one  needs  to  switch  from  the 
penalty  to  the  averaging  procedure.  Moreover  this  happened  only  at  very  few  time  steps. 


4  Uncertainty  Analysis  in  Supersonic  Flow  Problems 


It  is  well  known  that  the  steady  state  of  an  isentropic  flow  in  a  double-throated  nozzle  is 
not  unique.  In  fact,  the  steady  state  flow  can  be  either  completely  supersonic  or  completely 
subsonic  or  a  flow  containing  a  shock  wave  connecting  the  supersonic  branch  of  the  solution 
to  the  subsonic  branch,  the  location  of  the  shock  wave  depends  uniquely  on  the  initial 
condition. 


In  [51]  a  model  equation,  having  the  same  features,  was  considered  and  analyzed,  and 
in  [11]  the  model  has  been  generalized. 

In  many  applications  there  are  uncertainties  involved  in  the  initial  conditions  and  the 
question  arises:  what  can  be  said  about  the  shock  location  if  there  are  uncertainties  in  the 
initial  conditions.  We  note  that  while  randomness  enters  through  the  initial  conditions  in 
this  problem,  random  effects  can  generally  enter  into  practical  problems  through  boundary 
conditions,  initial  conditions,  the  domain  geometry,  missing  variables  and  fluid  properties 
:  etc.  Such  random 

new  methodologies  to  model  and  analyze  the  impact  of  such  uncertainties.  In  our  case 
we  are  interested  in  the  statistics  of  some  derived  quantities  (e.g.,  the  shock  position  of  a 
solution).  Such  are  often  hard  to  accurately  compute  from  the  first  few  moments  of  the 
solutions. 


In  a  preliminary  work,  generalized  polynomial  chaos  methods  were  implemented  to 
compute  the  probability  density  function  (PDF)  of  the  shock  location  for  the  cases  where 
the  initial  conditions  are  assumed  to  be  different  random  processes  (fields). 

Our  preliminary  conclusions  are: 


•  Polynomial  chaos  (PC)  expansion  modes  are  smooth  functions  of  the  spatial  vari¬ 
able  x,  although  the  individual  solution  realizations  are  discontinuous  in  the  spatial 
variable  x. 


•  The  solution  is  discontinuous  in  the  random  variable  space  at  a  fixed  point  x.  Filtering 
is  necessary  for  the  stability  of  the  scheme,  because  generalized  polynomial  chaos 


methods  are  spectral  representations  of  the  random  processes. 

•  When  the  variance  of  the  initial  condition  is  small,  the  probability  of  the  density 
function  (PDF)  of  the  shock  locations  is  computed  with  high  accuracy.  Otherwise, 
many  PC  expansion  terms  are  needed  to  produce  reasonable  results.  As  first  noted 
by  Chorin,  this  is  due  to  the  slow  convergence  of  PC  expansions  of  discontinuous 
functions  in  the  random  fields. 

•  The  biggest  absolute  eigenvalue  of  the  Jacobi  matrix  of  the  system  increases  quickly 
with  respect  to  the  number  of  PC  terms  used  in  the  expansion.  This  might  cause 
large  dissipation  for  some  numerical  schemes.  The  fast  increasing  size  of  the  system, 
when  using  more  PC  terms,  could  also  be  problematic  if  one  wants  to  solve  the 
system  with  a  high  order  numerical  scheme  using  characteristic  decomposition,  e.g., 
high  order  ENO  or  WENO. 

The  fact  that  the  coefficients  in  the  PC  expansion  are  smooth  is  surprising.  It  seems 
although  the  solution  contains  a  shock,  we  can  compute  it  by  embedding  it  in  a  random 
space  and  computing  smooth  solutions  in  that  space.  We  plan  to  study  whether  this  might 
be  a  general  procedure  in  computing  shock  waves. 


5  Richtmyer-Meshkov  Instabilities 

In  the  mixing  of  fuel  with  oxidants  in  SCRAM  JET  engine,  the  fuel  jets  are  under  impulsive 

the  fuel-oxidants  interface  and  the  breakup  of  fuels  into  finer  droplets.  Inertial  Confinement 
Fusion  Program  (ICF)  uses  high  energy  pulse  source  such  as  X-ray  and  laser  to  illuminate 
the  target  sphere  in  order  to  achieve  auto  fusion  ignition  on  the  National  Ignition  Facility. 
It  is  crucial  to  achieve  an  uniform  compression  of  the  target  sphere  as  possible  for  maximum 
efficiency  as  the  impulsively  accelerated  non-uniform  sphere  surface  by  a  shock  wave  causes 
a  non-uniform  pressure  profile  over  the  sphere. 

The  source  of  these  phenomenon  known  as  Richtmyer-Meshkov  Instability  is  related  to 
the  well  known  fluid  instability  studied  theoretically  by  Richtmyer  and  experimentally  by 
Meshkov.  The  Richtmyer-Meshkov  Instability  (RMI)  results  from  a  impulsively  accelerated 
interface  of  materials  with  different  densities  under  perturbation.  This  form  of  instability  is 
different  than  the  closely  related  fluid  instability  known  as  the  Rayleigh-Taylor  instability 
in  which  the  material  interface  is  under  constant  acceleration  force  such  as  gravity.  The 
vorticity  generated  by  the  cross  product  of  the  pressure  gradient  and  the  density  gradient 
deforms  the  and  amplifies  the  interface  perturbation  and  grows  in  time.  The  penetration  of 
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the  heavier  fluid  into  lighter  fluid  form  spike  and  bubble  vice  versa.  Growth  of  the  interface 
amplitude  and  secondary  shear  instability  promote  the  onset  of  turbulence  mixing  and 
enlarged  in  time.  The  RMI  is  encountered  in  a  variety  of  physical  contexts  such  as,  but 
not  limited  to,  those  described  above.  Reader  are  referred  to  extensive  literature  available. 

In  order  to  capture  the  shock-interface  interaction  and  the  fine  scale  structures  within 
the  turbulence  mixing  zone,  high  order  methods  are  highly  desirable.  Among  the  high 
orders  schemes  considered,  Spectral  methods  (Spectral)  and  Weighted  Essentially  Non- 
Oscillatory  finite  difference  schemes  (WENO)  are  considered  in  this  study.  High  Order 
compact  scheme  is  another  candidate  but  was  not  considered  in  this  study.  High  order,  in 
the  sense,  means  order  of  accuracy  is  at  least  greater  than  two.  In  this  study,  we  devise  the 
algorithms  based  on  the  methods  mentioned  above.  To  our  knowledge,  this  is  probably  the 
first  time  these  high  order  methods  are  implemented  for  the  study  of  the  RMI.  In  order  to 
evaluate  the  performance  of  the  devised  schemes,  only  single  mode  interface  perturbation 
will  be  included  in  this  study. 

From  the  point  of  view  of  the  numerical  calculation,  we  can  break  the  RMI  problem 
into  two  parts. 

First,  we  have  the  issue  of  reliably  calculating  the  motion  of  a  possibly  very  strong 
shock  wave,  and  second  we  have  the  issue  of  reliably  calculating  the  mix  that  ensues  after 
this  shock  wave  accelerates  the  interface.  It  is  in  this  second  area  of  calculating  the  ensuing 
mix  where  high  order  numerical  schemes  offer  unparalleled  efficiency.  This  efficiency  comes 
from  the  very  fundamental  fact  that  the  truncation  error  in  the  differentiation  operators 
can  be  made  much  smaller  by  increasing  the  order  of  the  scheme  than  by  increasing  the 
number  of  grid  points  .thereby ^  makmg_  jong  jime=  integration_computationally  much  less 
expensive  with  high  order  schemes. 

Secondly,  when  the  flow  variables  are  differentiable  and  contain  significant  structure, 
i.e.,  when  a  Fourier  representation  of  the  flow  variables  decays  slowly,  then  the  truncation 
error  is  again  highly  favorable  to  high  order  schemes. 

For  these  two  reasons,  if  one  has  structure,  with  or  without  shocks,  high  order  schemes 
are  orders  of  magnitude  computationally  more  efficient  than  low  order  schemes.  This  leads 
us  to  the  issue  of  high  order  schemes  in  the  presence  of  strong  shocks.  Of  course,  if  one  has 
only  shocks  and  no  structure,  then  there  is  no  reason  for  using  high  order  schemes.  So,  we 
assume  the  existence  of  structure  in  the  flow  variables. 

When  shocks  are  present,  there  has  been  a  great  deal  of  technology  developed  over 
the  last  few  decades  that  insures  that  one  can  obtain  a  high  order  accuracy  away  from 
the  shock  even  though  the  calculation  is,  as  it  must  be,  first  order  at  the  shock.  Here  we 
explore  two  such  high  order  schemes,  the  Weighted  Essentially  Non-oscillatory  (WENO) 
scheme  and  spectral  methods.  The  WENO  scheme  and  its  predecessor  the  Essentially 
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Non-oscillatory  (ENO)  schemes  were  designed  exactly  for  problems  such  as  RMI  and  are 
extremely  robust  even  for  very  strong  shocks.  We  would  like  to  stress  the  importance  of 
comparing  numerical  results  obtained  by  algorithms  based  on  different  philosophy.  The 
confidence  of  the  numerically  converged  results  can  be  greatly  enhanced  by  their  agreements 
of  the  numerical  results.  In  our  case,  spectral  methods  is  global  in  nature  as  opposed  to 
the  finite  difference  scheme  which  is  local  in  nature  as  seen  in  the  following  sections. 

In  a  series  of  papers  [32,  33,  34,  35,  36],  it  was  shown  that  given  the  spectral  approx¬ 
imation  to  a  piecewise  smooth  function  one  can  construct  an  exponentially  convergent 
approximation  to  the  function.  This  result  had  been  proven  for  Fourier,  Chebyshev,  Leg¬ 
endre  and  spectral  methods  based  on  spherical  harmonics. 


For  the  ease  of  presentation,  we  discuss  the  results  for  Fourier  approximations,  since 
nothing  essential  is  lost  in  the  Chebyshev  case. 

The  exponential  filter  offers  the  flexibility  of  changing  the  order  of  the  filter  simply  by 
specifying  a  different  7.  One  does  not  have  to  write  a  different  filter  for  different  order. 
Thus  varying  7  with  N  yields  exponential  accuracy  according  to  [59]. 

The  above  defined  filters  do  not  completely  remove  the  Gibbs  phenomenon  as  oscil¬ 
lations  still  exist  in  the  neighborhood  of  the  discontinuities.  In  order  to  recover  the  full 
accuracy  in  any  region  where  the  function  is  continuous,  one  has  to  use  a  different  idea.  In 
[35]  it  is  shown  how  to  use  a  known  set  of  2N  +  1  Fourier  coefficients  to  obtains  the  coeffi¬ 
cients  of  a  different  expansion  (based  on  the  Gegenbauer  polynomials).  The  new  expansion 
converges  exponentially  in  any  smooth  region. 


In  practice,  when  solving  differential  equations  one  uses  the  exponential  filter  at  every 

time  step  and  the -Gegenbauer  filter  at  the  end  of  the  calculations  as^  a  postprocessor.  *  ~ 


When  spectral  methods  are  applied  to  nonlinear  hyperbolic  equations  in  conservation 
form  the  problem  of  an  entropy  satisfying  solution  arises.  In  fact,  there  is  no  artificial 
dissipation  in  the  method  to  indicate  that  the  solution  is  a  limit  of  a  dissipative  process. 
One  clearly  needs  to  add  artificial  dissipation.  However  this  dissipation  must  be  spectrally 
small  in  order  not  to  affect  the  overall  accuracy.  This  problem  had  been  addressed  in 
[56,  57,  44], 

It  has  been  shown  that  with  a  suitable  addition  of  (spectrally  small)  artificial  dissipation 
to  the  high  modes  only,  the  method  converges.  In  this  paper  we  used  one  version  of  the 
above  idea:  the  Super  Vanishing  Viscosity  method  (SVV)  suggested  by  Tadmor  [56,  57]. 

The  interaction  between  the  shock  and  the  interface,  generates  a  shock  triple  point 
along  the  gaseous  interface.  Localized  sharp  gradient  could  cause  numerical  instability  if 
insufficient  localized  physical/numerical  dissipation  existed  there.  For  the  Spectral  scheme, 
instabilities  were  observed  in  the  case  of  under- resolved  simulation  and/or  using  too  high 
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order  of  the  global  filtering.  Global  strong  filtering  would  keep  the  scheme  stable  at  the 
cost  of  smoothing  out  all  fine  scale  physical  structures,  which  is  highly  undesirable  when 
resolution  of  fine  scale  structures  is  essential  for  understanding  of  the  issue  of  shock  induced 
mixing. 

Since  this  is  a  local  phenomenon,  a  sufficient  strong  local  dissipation  or  Local  Adaptive 
Filtering  would  be  needed  to  keep  the  scheme  stable.  The  location  of  those  collocation 
points  where  the  density  p  <=  pt0i  will  first  be  located  and  marked  for  further  processing. 
A  global  strong  filtering,  for  example  7  =  2  will  be  used  to  reduce  the  magnitude  of  the 
oscillations  at  those  points  only.  In  all  the  cases  studied  here,  ptoi  is  set  to  be  \pmin  where 
pmin  is  the  expected  minimal  density  value,  pm,n  =  Pxe  for  example. 

Our  experiences  with  this  class  of  problem  indicated  that  the  local  adaptive  filtering 
only  applied  on  a  few  (in  the  range  of  1  to  7)  grid  points  around  the  shock  triple  point 
only  as  the  shock  propagated  along  the  gaseous  interface  in  the  earlier  time.  Otherwise, 
the  Spectral  scheme  remained  stable  for  all  other  time  and  no  local  adaptive  filtering  is 
needed. 

After  the  evaluation  of  the  spatial  operator,  the  third  order  TVD  Runge-Kutta  scheme 
of  Shu  and  Osher  [54]  was  used  to  solve  the  system  of  ODE’s.  It  has  the  form  of 


U 1  =  Un  +  AtL(Un) 

U 2  =  i(3t/n  +  C/1  +  AtL{U1)) 

Un+1  =  ~(Un  +  2U2  +  2AtL(U2)). 

O 


L  is  the  spatial  operator.  UnmS  Un+x~sxe  the ‘data  arrays  at  the  rr  th  and  (h  f  lj-th  time 
step,  respectively.  L^and  U2  are  two  temporary  arrays  at  the  intermediate  Runge-Kutta 
stages.  The  scheme  is  stable  for  CFL  <  1.  This  Runge  Kutta  scheme  is  a  low  storage 
scheme  as  it  can  be  rewritten  in  such  a  way  that  only  two  arrays  Un  and  Ul  are  needed. 


We  summarize  the  overall  solution  procedures  of  the  algorithm  as  follow, 


•  Periodical  Domain  is  specified  in  the  y  direction  and  symmetry  property  of  the  prob¬ 
lem  is  exploited  to  reduce  the  amount  of  computational  operations  by  half. 

•  Spatial  Algorithm  : 

1.  Combined  Chebyshev  ( x )  and  Fourier  ( y )  collocation  method  (Spectral), 

—  Differentiation  and  Smoothing  operations  are  performed  via  an  optimized 
library  PseudoPack  [13]; 

-  A  10’th  and  9’th  order  exponential  filter  used  for  the  differentiation  and 
solution  smoothing  respectively; 


—  Kosloff-Tal-Ezer  mapping  [24,  21]  for  the  Chebyshev  collocation  methods. 

-  An  adaptive  local  smoothing. 

2.  Fifth  order  WENO  finite  difference  scheme  (WENO-LF-5)  with  Lax-Frederick 
flux  splitting. 

•  The  third  order  TVD  Runge  Kutta  method  by  Shu  and  Osher  [54]  is  used  to  advance 
the  solution  in  time. 

A  series  of  numerical  simulations  are  carried  out  to  investigate  the  convergence  prop¬ 
erties  of  both  the  Spectral  scheme  and  the  WENO  scheme.  Simulations,  using  various 
interface  thicknesses  and  resolutions,  are  computed  and  terminated  at  some  representative 
time  after  the  shock  had  transmitted  sufficiently  far  away  from  the  interface  and  before 
exiting  the  physical  domain.  It  allows  the  development  of  vortical  rollups  of  the  gaseous 
interface.  Vorticities  are  generated  by  the  cross  product  of  the  pressure  gradient  of  the 
shock  and  the  density  gradient  of  the  gases.  The  final  time  is  set  to  t  =  50  x  10~6s  for 
Lx  =  50  cm  and  t  =  143  x  10_6s  for  Lx  =  150  cm. 

As  evidenced  from  the  results  of  the  Spectral  and  the  WENO  calculations  shown  below, 
the  following  major  features  of  the  Richtmyer-Meshkov  instability  can  be  observed  (see 
figure  3)  at  time  t  =  50  x  10~6s,  namely, 

•  Wave  generated  by  the  shock  refraction  behind  the  gas  interface  in  Box  1. 

•  The  penetration  of  the  heavy  (Xe)  to  light  (Ar)  fluid  causes  the  deformation  of  the 
interface  into ■?;  large  ipushroom  shape  structures-in Box  2  .and  the.opposite  in  Box  5. 
They  are  referred  as  Spike  and  Bubble  respectively,  in  the  literatures.  They  move  in 
the  opposite  direction  relative  to  each  other  and  form  a  ever  larger  turbulence  mixing 
zone. 

•  Pressure  wave  along  the  transmitted  shock  in  Box  3. 

•  A  small  jet  and  its  vortical  structure  located  in  Box  4.  The  contact  discontinuity 
develops  into  a  more  complicated  vortical  rollups  in  a  finer  and  long  term  simulation 
possibly  caused  by  the  Kelvin-Helmholtz  instability. 

•  Vortical  rollups  of  the  gaseous  interface  inside  Box  6. 

The  global  large  and  median  features  (Box  1,  2,  3,  4  and  5)  are  well  captured  accurately 
by  both  numerical  schemes  for  a  given  resolution.  It  is  unclear,  however,  if  the  smaller 
rollups  along  the  gases  interface  (Box  6)  presented  in  the  high  resolution/high  order  cases 
are  physical  due  to  the  non-dissipative  nature  of  the  Euler  equations  or  numerical  due  to 
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Figure  1:  The  numbered  regions  enclose  the  most  prominent  flow  features  of  the  Richtmyer- 
Meshkov  instability  at  time  t  —  50  x  1CT6s. 

the  oscillatory  nature  of  the  numerical  schemes  or  both.  More  researches  are  needed  to 
answer  this  question  fully. 

We  shall  first  examine  the  convergence  property  of  both  the  Spectral  scheme  and  the 
WENO-LF-5  finite  difference  scheme.  For  this,  we  used  a  thicker  interface  with  S  = 
0.6  cm  to  establish  the  convergence  of  the  numerical  schemes  of  the  large  and  medium 
scale  structures  (box  1,  2,  3,  4  and  5  in  figure  3).  This  avoids  the  possible  contamination 
of  numerical  artifacts  due  to  high  gradients  generated  along  the  shock-interface  interaction 
and  bypass  the  issue  of  under-resolved  fine  scale  physical  structures.  Furthermore,  the 
spectral  solutions  are  not  post-processed  by  any  existing  post-processing  algorithms  to 
remove  the  Gibbs  oscillations. 

Convergence  Stu^y:  :  6  =JL6 cm  _  4  . . . 

The  density  p  and  velocity  V  of  the  solution  of  the  Spectral  and  WENO-LF-5  runs  are 
shown  in  figure  4  at  time  t  =  50  x  10-6  s  with  various  resolutions. 

It  can  observed  that  the  large  and  medium  scale  structures  such  as  the  transmitted 
shock,  the  location  of  the  triple  point,  the  shocked-interface  velocity,  pressure  waves  and 
vorticity  generation,  are  basically  in  excellent  agreement  with  each  others.  The  weak 
vertical  wave  located  downstream  behind  the  interface  is  an  left  over  entropy  wave  from 
the  initial  shock  condition. 

Convergence  Study  :  5  =  0.2  cm 

The  interface  thickness  is  further  reduced  significantly  from  6  =  0.6  cm  to  5  =  0.2  cm. 
The  density  p  and  velocity  V  of  the  solution  of  the  Spectral  and  WENO-LF-5  runs  are 
shown  in  figures  5  and  6  respectively,  at  time  t  =  50  x  10-6  s  with  various  resolutions. 

Similar  to  the  previous  case  of  5  —  0.6  cm,  it  can  be  observed  that  the  large  and  median 


Figure  2:  Convergence  Study  5  =  0.6  cm  :  Density  (Top  Row)  and  V- Velocity  (Bottom 
Row)  contour  plot  of  the  Richtmyer- Meshkov  instability  as  computed  by  the  Spectral 
scheme  and  the  WENO-LF-5  scheme.  Domain  length  in  x  is  Lx  =  5  cm.  The  interface 
thickness  6  —  0.6  cm.  The  final  time  is  t  =  50  x  10~6s.  The  resolution  of  the  Spectral 
schemes  are  256x128  (Top  Left),  512x256  (Top  Right)  and  1024x512  (Bottom  Left)  and 
the  WENO  scheme  is  1024x512  (Bottom  Right). 


scale  structures  such  as  transmitted  shock,  shocked-interface  velocity  and  shock  triple  point 
are  basically  in  excellent  agreement  with  each  others.  Some  discrepancies  of  the  fine  scale 
structures  along  the  gaseous  interface,  as  can  be  expected  for  numerical  simulation  of  the 
Euler  equations . sensitive , to  perturbation  in  naturerare  observed.  ... . 

Snapshot  of  the  evolution  of  density  and  velocity  flow  fields  at  several  immediate  times 
are  illustrated  in  figure  (7),  for  the  Spectral  scheme  and  in  figure  (7)  for  the  WENO  scheme. 
The  contour  levels  are  the  same  and  constant  for  both  schemes  in  all  plots. 

The  Mach  number  M,  the  Atwood  number  At  and  the  interface  curvature  play  an 
important  role  on  the  growth  of  perturbed  amplitude  on  the  interface.  In  the  particular 
set  of  parameters  studied  here  with  high  Mach  number  M  =  4.46  and  median  Atwood 
number  At  ~  0.54,  a  formation  of  triple-shock  configuration  along  the  interface  indicates 
that  shock-interface  interaction  is  in  the  "hard”  regime.  A  ’’hard”  regime,  as  quoted  from 
Zaytsev  etc.  is  ’’the  propagation  of  secondary  shocks  across  the  flow  that  is  accompanied 
by  the  formation  of  breaks  and  triple  configurations  on  the  refracted  and  reflected  shocks” . 
The  triple-shock  formation  can  be  observed  easily  in  the  early  time  <~  30  x  10_6s. 

As  the  Mach  number  increases,  high  order  schemes  tend  to  break  down  as  the  solution 
develops  singularity  and/or  high  gradients  in  time  in  which  the  numerical  schemes  becomes 
more  difficult  to  resolve  them.  For  the  spectral  methods,  we  have  developed  a  Spectral 
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Figure  3:  Convergence  Study  5  =  0.2  cm  :  Density  (Top  Row)  and  V- Velocity  (Bottom 
Row)  contour  plot  of  the  Richtmyer-Meshkov  instability  as  computed  by  the  Spectral 
scheme.  Domain  length  in  x  is  Lx  =  5  cm.  The  interface  thickness  8  =  0.2  cm.  The  final 
time  is  t  =  50  x  10“6s.  The  resolution  of  the  Spectral  schemes  are  384x192  (Top  Left), 
512x256  (Top  Right)  and  1024x256  (Bottom  Left). 

Adaptive  Domain  Algorithm  which  adjusts  the  size  of  the  computational  domain  in  time 
in  order  to  resolve  the  fine  scale  structures  at  high  Mach  number  as  they  are  developing  in 
time.  In  figure  (8),  we  presents  the  density  for  the  Richtmyer-Meshkov  instability  for  the 
case  similar  to  the  previous  section  except  that  the  Mach  number  is  increased  from  Mach 
4.46  to  Mach -10V, ?*:*■_  :-y.  ;r: r  :sr 

In  figure  9  we  show  the  preliminary  result  of  the  density  isosurface  plot  of  the  RMI 
with  a  Mach  4.46  shock  and  a  three  dimensional  random  pertubation  interface  separating 
the  Argon  and  Xenon  gases. 


6  Recessed  Cavity  Flameholders 

Recessed  cavities  provide  a  high  temperature,  low  speed  recirculating  region  that  can  sup¬ 
port  the  production  of  radicals  created  during  chemical  reactions.  This  stable  and  efficient 
flame-holding  performance  by  the  cavity  is  achieved  by  generating  a  recirculation  region 
inside  the  cavity  where  a  hot  pool  of  radicals  forms  resulting  in  reducing  the  induction  time 
and  thus  obtaining  the  auto-ignition  [6,  61].  Experiments  have  shown  that  such  efficiency 
depends  on  the  geometry  of  the  cavity  such  as  the  degree  of  the  slantness  of  the  aft  wall  and 
the  length  to  depth  ratio  of  cavity  L/D.  Thus  one  can  optimize  the  flame-holding  perfor¬ 
mance  by  properly  adjusting  the  geometrical  parameters  of  the  cavity  flame-holder  system 
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Figure  4:  Convergence  Study  6  =  0.2  cm  :  Density  (Top  Row)  and  V-Velocity  (Bottom 
Row)  contour  plot  of  the  Richtmyer- Meshkov  instability  as  computed  by  the  WENO-LF- 
5scheme.  Domain  length  in  x  is  Lx  =  5  cm.  The  interface  thickness  6  =  0.2  cm.  The 
final  time  is  t  =  50  x  10_6s.  The  resolution  of  the  WENO-LF-5  schemes  are  256x128  (Top 
Left),  512x256  (Top  Right)  and  1024x512  (Bottom  Left). 


for  a  given  supersonic  flight  regime.  There  are  two  major  issues  of  such  cavity  flame-holder 
system  that  need  to  be  investigated  ;  (1)  What  is  the  optimal  angle  of  the  aft  wall  for  a 
given  L/D ?  and  (2 )How  does  the  fuel  injection  interact  with  cavity  flows?  An  answer  to 
these  questions  require  both  a  comprehensive  laboratory  and  numerical  experiments. 


-  Results  of  sS^^raf  SufficncM  sCuHfes^Have  shown  that"  thest  ability  of  the  recirculation 
inside  cavity  is  enhanced  for  the  lower  angle  of  cavity  compared  to  the  rectangular  cavity. 
The  present  study,  however,  gives  more  accurate  and  finer  details  of  the  fields  than  those 
done  by  lower  order  numerical  experiments.  We  show  that  a  stationary  recirculation  region 
is  not  formed  inside  the  cavity  contrary  to  what  the  lower  order  schemes  predict.  A 
quantitative  analysis  made  in  this  study  shows  that  the  lower  angled  wall  of  the  cavity 
reduces  the  pressure  fluctuations  significantly  inside  the  cavity  for  the  non-reactive  flows. 
We  obtained  a  similar  result  for  the  reactive  flows  with  the  ignition  of  the  fuel  supplied 
initially  in  the  cavity. 


In  the  SCRAM  Jet  community,  a  cavity  with  the  length-to-depth  ratio  L/D  <  7  ~  10 
is  usually  categorized  as  an  ’open’  cavity  since  the  upper  shear  layer  re-attaches  at  the 
back  face  [6].  In  this  work,  we  choose  the  L/D  of  the  baseline  cavity  to  be  4  and  thus  the 
open  cavity  system  is  considered.  The  coordinates  of  the  cavity  are  (7 cm,  —1cm)  for  the 
upper  left  and  (11cm,  —2 cm)  for  the  right  bottom  corners  of  cavity.  With  the  length  of  the 
neck  of  the  cavity  fixed  to  be  4cm,  we  consider  three  different  angles  of  the  right  corner 
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Figure  5:  Snapshot  of  evolution  with  the  Spectral  (Leftmost  two  figures)  and  WENO-LF-5 
(Rightmost  two  figures)  schemes  and  S  =  0.2  cm  : 

In  the  leftmost  two  figures  (Spectral),  the  Density  (Left)  and  Velocity  (Right)  contour 
plot  of  the  Richtmyer-Meshkov  instability  as  computed  by  the  Spectral  scheme  at  time 
t  =  12.5  x  10~6s,  25.0  x  10-6s,  31.3  x  10-6s,  37.5  x  10"6s,  43.8  x  10~6s  and  50.0  x  10"6s  are 
showed.  The  resolution  of  the  Spectral  scheme  is  1024x256. 

In  the  rightmost  two  figures  (WENO-LF-5),  the  Density  (Left)  and  Velocity  (Right)  con¬ 
tour  plot  of  the  Richtmyer-Meshkov  instability  as  computed  by  the  WENO-LF-5  scheme  at 
time  t  =  13.0  x  10"6s,  24.7  x  10~6s,  31.5  x  10~6s,  37.1  x  10"6s,  43.2  x  10“6s  and  50.0  x  10"6s 
are  shown.  .The  r  esoJution  of  the  WENO-LF-5  scheme  is!024x512.  -  --  — 

Domain  length  in  x  is  Lx  =  5  cm.  The  interface  thickness  5  =  0.2  cm. 

of  the  floor  of  the  cavity  (  60,45  and  30),  we  then  compare  each  one  with  the  case  of  the 
rectangular  aft  wall.  The  fluid  conditions  are  given  as  followings;  the  free  stream  Mach 
number  M  =  1.91,  total  pressure  P  —  2.82 (atm),  total  temperature  T  =  830.6(iQ  and 
normalized  Reynolds  number  Re  =  3.9  x  107(l/m).  Note  that  the  Reynolds  number  is  here 
normalized  and  has  a  unit  of  1  /[length],  also  the  Reynolds  number  based  on  the  cavity 
dimensions  is  O(105).  The  boundary  layer  thickness  scale  is  6  =  5  x  10_4(m),  and  finally, 
the  wall  temperature  is  Tw  =  460.7835(A').  The  initial  configuration  for  the  baseline  cavity 
system  is  shown  in  figure  10. 

We  have  conducted  two  different  experiments  for  each  of  the  following  cases  :  (1)  non¬ 
reacting  cold  flow ,  and  (2)  reacting  flow  .  We  use  9  and  17  subdomains  for  both  cases  1 
and  2.  For  the  outflow  conditions  at  the  exit  of  the  system  and  at  the  upper  boundary, 
we  mainly  use  a  semi-infinite  mapping  in  order  to  reduce  the  possible  reflections  at  the 
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Figure  6:  High  Mach  number  case  (5  =  0.2  cm  :  Density  contour  plot  of  the  Richtmyer- 
Meshkov  instability  as  computed  by  the  Spectral  scheme.  The  Mach  number  M=10.  Do¬ 
main  length  in  x  is  Lx  =  30  cm.  The  interface  thickness  8  =  0.2  cm.  The  final  time  is 
t  —  125  x  10_6s.  The  resolution  of  the  Spectral  schemes  are  2304x256. 

boundaries.  The  characteristic  boundary  conditions  are  also  applied  and  will  be  discussed 
in  the  next  section  and  compared  to  the  mapping.  For  the  case  of  the  reactive  flows,  the 
cavity  was  initially  filled  with  Hydrogen  fuel  with  fuel-to-total  gas  ratio  of  0.5.  The  order  of 
the  polynomial  of  approximation  in  y  direction  in  the  domain  beside  the  wall  is  taken  large 
enough  to  resolve  the  boundary  layer  well.  Finally  the  adaptive  filtering  is  turned  on  if  the 
mass  fraction  of  Hydrogen  and  Oxygen' exceed  the  range  of  —  0.09  <  fn2  £  1.09,  —0.02  <  ' 
fo2  <  0.25  and  the  temperature  exceeds  the  range  of  300(JT)  <  T  <  3500(JT).  As  the  shear 
layer  and  the  complex  features  of  the  flows  develop,  the  adaptivity  criteria  for  applying 
the  local  smoothing  is  satisfied  at  some  points.  In  the  calculations,  we  use  the  3rd  and  2nd 
order  local  filtering  for  the  non-reactive  and  reactive  flows  respectively.  It  turns  out  that 
the  local  smoothing  was  applied  in  very  few  points  at  the  upper  corner  of  the  cavity  wall. 

For  the  adaptive  averaging,  we  use  the  criteria  constant  Cave  such  that  the  difference 
of  the  state  vectors  (or  pressure)  between  the  two  adjacent  domains  is  less  than  10%.  In 
figure  11  the  Penalty  Navier-Stokes  equations  were  considered  for  the  non-reactive  cold 
flows.  As  evident  from  the  contours  of  the  density,  the  approximations  were  well  matched 
at  the  interfaces.  Here  the  outer  boundary  was  approximated  by  using  the  characteris¬ 
tic  conditions  of  the  inviscid  fluxes.  The  adaptive  averaging,  with  the  given  adaptivity 
conditions  above,  took  place  at  only  a  few  points.  The  characteristic  boundary  condi¬ 
tions  using  the  inviscid  fluxes  yield  good  results  for  both  the  problems  of  the  density  peak 
propagation  and  the  non- reactive  cold  flows.  As  in  figure  1,  we  observe  that  there  exist 


Figure  7:  Density  isosurface  plot  of  the  RMI  with  a  Mach  4.46  shock  and  a  three  dimen¬ 
sional  random  perturbation  interface.  The  gases  are  the  Argon  and  Xenon 
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Figure  8:  The  initial  configuration  for  the  baseline  cavity  system. 

penalty  parameters  satisfying  the  stability  conditions  that  may  induce  reflecting  modes  at 
the  interfaces.  Figure  12  shows  the  pressure  history  of  the  non-reactive  cold  flows  for  the 
various  angles  of  the  aft  wall  at  two  different  locations  inside  the  cavity,  i.e.  at  the  center, 
(x,y)  =  (8.5 cm,  —1.5 cm),  and  at  the  middle  of  the  floor  (x,  y)  =  (8.5 cm,  —1.9 cm). 

These  figures  show  that  the  pressure  fluctuations  in  cavities  with  lower  angle  of  the  aft 
are  weaker  than  in  cavities  with  higher  angles.  It  is  also  shown  that  the  attenuation  of 
the  pressure  fluctuations  are  obtained  both  at  the  center  and  the  middle  of  the  floor  of 
the  cavity.  It  is  interesting  to  observe  that  the  patterns  of  the  pressure  fluctuations  for 
a  given  angle  at  different  locations  are  different  depending  on  the  angle.  In  the  case  of 
the  30  degree  aft  wall,  the  pressure  fluctuations  are  almost  the  same  at  the  two  locations 
considered  whereas  the  case  of  45  degree  shows  a  difference  in  the  patterns  of  the  pressure 


Figure  9:  The  non-reactive  cold  flows  with  the  penalty  Navier-Stokes  equations:  the  density 
contours  are  given  in  this  figure  at  t  =  0.25ms.  17  domains  are  used  and  the  boundaries 
of  each  domain  are  shown. 
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Figure  10:  Pressure  history  for  non-reactive  flows:  the  left  panel  represents  the  pressure 
history  at  the  center  of  the  cavity  and  the  right  panel  at  the  middle  of  the  floor  of  the  cavity. 
Each  panel  shows  the  case  of  90,  60,  and  30  degree  cavity  walls  from  top  to  bottom. 

fluctuations  between.  the^twq.Jocations,.  .Jhe  ^pressure  fluctuations  a^  the  bottom  grows 
greater  than  that  at  the  center  after  some  time. 

Figure  13  shows  the  pressure  history  when  the  heavy  global  filter  is  applied  (in  this  case, 
the  4th  order  filter  was  used).  Unlike  the  previous  case  illustrated  in  figure  12,  where  the 
6th  order  global  filter  is  used,  the  pressure  fluctuations  eventually  decay  out  and  a  large 
recirculation  zone  is  formed  inside  the  cavity  without  any  severe  pressure  fluctuations. 
Note  that  the  scale  in  the  left  panel  shown  is  the  same  as  in  figure  12  while  the  right  panel 
is  shown  in  a  smaller  scale  for  a  closer  look.  This  figure  shows  that  the  large  recirculation 
zone(s)  formed  inside  the  cavity  obtained  by  the  lower  order  numerical  scheme  is  induced 
not  physically  but  rather  artificially  due  to  the  heavy  numerical  dissipations.  This  is  clearly 
shown  in  figure  14.  In  this  figure  a  large  recirculation  zone  is  observed  -  this  zone  is  formed 
earlier  than  this  streamlines  are  captured  -  when  the  4th  order  filter  is  used(left  figure) 
and  an  almost  steady  state  is  already  reached  as  the  pressure  history  indicates  in  figure 
13.  We  find  from  the  numerical  results  that  the  large  recirculation  is  very  stable  once  it 
forms.  This  large  recirculation  and  the  steady  state  solutions  are  not  observed  in  the  case 
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Figure  11:  Pressure  history  of  the  non-reactive  flows  with  the  use  of  the  fth  order  filter:  the 
left  panel  represents  the  pressure  history  at  the  center  of  cavity  and  the  right  panel  shows 
the  ■ left  panel  in  a  smaller  scale.  Each  panel  shows  the  case  of  90,  and  30  degree  cavity 
walls  from  top  to  bottom.  Note  that  the  scale  of  the  right  panel  is  different  from  the  left. 


Figure  12:  Streamlines:  the  left  figure  shows  the  streamlines  at  t  =  1.685ms  for  the  global 
filtering  order  7  =  4  and  the  right  at  t  =  2.38 ms  for  7  =  6. 

of  7  =  6(right).  For  the  case  of  7  =  6  instead  of  the  large  single  recirculation  zone,  smaller 
scale  vortex  circulations  are  formed  and  they  are  interacting  with  each  other,  never  reaching 
the  steady  state  with  time.  This  result  shows  that  for  these  sensitive  problems,  high  order 
accuracy  should  be  used  in  order  to  minimize  the  effect  of  the  numerical  dissipation. 

Figure  15  shows  the  case  of  the  reactive  flows  for  the  90  and  30  degree  aft  walls.  Similar 
•features  of- the  pressure-fluctuations 'are  shown-as  in  themon-reactiye  flows.-  However  the - 
pressure  fluctuations  are  much  more  attenuated  for  both  the  90  and  30  degree  walls  than 
in  the  non-reactive  cold  flows.  In  the  reactive  cases  Hydrogen  fuel,  which  was  initially 
supplied  inside  the  cavity  was  consumed.  As  time  elapses,  the  fuel  is  consumed  out  with 
the  production  of  the  water  for  these  cases. 
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Figure  13:  Pressure  history  for  reactive  flows:  the  left  panel  represents  the  pressure  history 
at  the  center  of  cavity  and  the  right  panel  at  the  middle  of  the  floor  of  cavity.  Each  panel 
shows  the  case  of  90  and  30  degree  cavity  walls  from  top  to  bottom. 

These  results  demonstrate  that  simulations  of  cold  flows  do  not  necessarily  shed  light 
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Figure  14:  The  water  contour  of  the  reactive  flows:  the  left  the  water  density  contour  is 
given  in  the  left  figure  and  its  streamlines  in  the  right  figure  at  t  =  0.135ms. 

on  the  behavior  of  reactive  flows. 

Figure  17  shows  the  density  contours  and  streamlines  for  the  90,  60,  45  and  30  degree 
walls  at  the  instant  time  t  =  2.4 ms.  As  shown  in  the  figure,  the  shear  layer  is  becoming 
weaker  as  the  degree  of  angle  of  the  aft  wall  and  the  flow  fields  are  becoming  more  regular¬ 
ized  for  the  case  of  the  lower  angle.  And  note  that  the  density  compression  at  the  corner 
of  the  aft  wall  is  also  becoming  weaker  for  the  more  slanted  wall  cases. 

Reactive  flow 

Figure  18  shows  the  water  contour  inside  the  cavity  for  the  different  angles  at  different 
time.  Here  we  define  the  region  where  the  flames  are  generated  to  be  same  as  the  region 
where  the  water  is  produced.  As  the  Hydrogen  fuel  is  consumed,  the  water  is  produced  and 
starts  to  be  expelled  from  the  cavity  to  the  main  channel.  The  flame-holding  efficiency  is 
enhanced  if  the  chemical  radicals  (water  in  this  case)  are  stably  circulating  and  long  lasting 
before  they  are  expelled  from  the  cavity.  Figure  18  shows  that  the  lower  angled  aft  wall 
(30  degree  in  this  case)  maintains  more  water  than  the  90  degree  wall  at  a  given  time.  The 
,  figure  also  shows  that,  the  lower  angled  aft  wall  holds  the  flame,  (water  Jn  this  ^ascj  jpnger  , 
than  the  90  degree  wall  -  in  the  last  figure  in  figure  18  at  t  =  2.26ms,  the  most  water  is 
expelled  and  the  only  the  small  amount  is  left  in  the  left  corner  while  the  30  degree  wall 
cavity  holds  the  water  still  through  the  cavity.  These  results  imply  that  the  flame-holding 
efficiency  can  be  increased  by  lowering  the  angle  of  the  aft  wall  of  the  cavity. 

Figure  19  shows  the  streamlines  corresponding  to  the  each  case  of  figure  18.  Note  that 
compared  to  the  non-reactive  cases,  the  shear  layers  are  less  developed  for  the  reactive 
cases.  As  the  figures  of  the  pressure  fluctuation  history  and  figure  19  indicate,  the  shear 
layers  are  weak  for  both  the  90  and  the  30  degree  walls  in  the  reactive  cases. 
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Figure  15:  The  density  contour  and  the  streamline  of  the  non-reactive  flows:  the  left  column 
shows  the  density  contour  for  90,  60,  45  and  30  degree  walls  from  top  to  bottom  and  the 
r  right  column  shows  the  corresponding -streamlines  at  t  =  2.43ms,  The  maximum  contour 
level  is  1.8  and  the  minimum  0.5  with  the  level  step  size  50. 
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